Learning vs Memory

Was the response learning or memory?

mod_E1 <-  rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
                                                                                        ~1|ES_ID,
                                                                                        ~1|Strain),
                  test = "t",
                  data = dat)

summary(mod_E1) 

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -9.3793   18.7586   28.7586   41.2577   29.4729   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0076  0.0873     30     no  Study_ID 
## sigma^2.2  0.0340  0.1843     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 90) = 802.5794, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 18.2861, p-val < .0001
## 
## Model Results:
## 
##                             estimate      se    tval  df    pval   ci.lb 
## Learning_vs_memoryLearning    0.2303  0.0450  5.1227  90  <.0001  0.1410 
## Learning_vs_memoryMemory      0.1624  0.0355  4.5713  90  <.0001  0.0919 
##                              ci.ub 
## Learning_vs_memoryLearning  0.3197  *** 
## Learning_vs_memoryMemory    0.2330  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_E1) 

##   R2_marginal R2_coditional 
##    0.02572583    1.00000000

LvsM_E<- orchard_plot(mod_E1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  xlim(-0.5, 2) + 
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21) + # mean estimate
  scale_size_continuous(range = c(1, 7)) + # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

LvsM_E

Fig. 4a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_E1 <- contrast_fun(data = dat, response = lnRR_Ea, moderator = Learning_vs_memory,
    VCV = VCV_E)
res_table_mod_E1 <- get_estimate(model = mod_E1, contra = contra_mod_E1, moderator = Learning_vs_memory)

res_table_mod_E1

## # A tibble: 3 x 7
##   Levels          Estimate Lower_CI Upper_CI    P_value Lower_PI Upper_PI
##   <chr>              <dbl>    <dbl>    <dbl>      <dbl>    <dbl>    <dbl>
## 1 Learning          0.230    0.141    0.320  0.00000170   -0.185    0.645
## 2 Memory            0.162    0.0919   0.233  0.0000154    -0.249    0.574
## 3 Learning-Memory  -0.0679  -0.164    0.0284 0.165        -0.484    0.349

Type of assay

The broad category of the type of assay used to measure learning or memory

dat1 <- filter(dat, Type_assay %in% c("Recognition", "Habituation", "Conditioning"))
VCV_E1 <- impute_covariance_matrix(vi = dat1$lnRRV_E, cluster = dat$Study_ID, r = 0.5)

mod_E2 <- rma.mv(yi = lnRR_Ea, V = VCV_E1, mod = ~Type_assay-1, random = list(~1|Study_ID,
                                                                                  ~1|ES_ID,
                                                                                  ~1|Strain),
                 test = "t",
                 data = dat1)

summary(mod_E2)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -7.8496   15.6993   27.6993   42.6311   28.7237   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0099  0.0995     30     no  Study_ID 
## sigma^2.2  0.0321  0.1792     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 661.8903, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.7993, p-val < .0001
## 
## Model Results:
## 
##                         estimate      se    tval  df    pval    ci.lb   ci.ub 
## Type_assayConditioning    0.2167  0.0351  6.1650  89  <.0001   0.1468  0.2865 
## Type_assayHabituation     0.1821  0.1147  1.5886  89  0.1157  -0.0457  0.4100 
## Type_assayRecognition     0.0554  0.0640  0.8659  89  0.3889  -0.0717  0.1826 
##  
## Type_assayConditioning  *** 
## Type_assayHabituation 
## Type_assayRecognition 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_E2) 

##   R2_marginal R2_coditional 
##    0.07376134    1.00000000

Learning_E <- orchard_plot(mod_E2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  xlim(-0.5, 2) + 
  scale_colour_manual(values = c("grey34","grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34", "grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
    geom_signif(comparisons = list(c("Conditioning", "Recognition")),
              annotations = "*") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Learning_E

Fig. 4b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_E2 <- contrast_fun(data = dat1, response = lnRR_Ea, moderator = Type_assay,
    VCV = VCV_E1)
res_table_mod_E2 <- get_estimate(model = mod_E2, contra = contra_mod_E2, moderator = Type_assay)

res_table_mod_E2

## # A tibble: 6 x 7
##   Levels                Estimate Lower_CI Upper_CI     P_value Lower_PI Upper_PI
##   <chr>                    <dbl>    <dbl>    <dbl>       <dbl>    <dbl>    <dbl>
## 1 Conditioning            0.217    0.147    0.287      2.02e-8   -0.197    0.630
## 2 Habituation             0.182   -0.0457   0.410      1.16e-1   -0.285    0.649
## 3 Recognition             0.0554  -0.0717   0.183      3.89e-1   -0.371    0.482
## 4 Conditioning-Habitua…  -0.0345  -0.262    0.193      7.63e-1   -0.501    0.432
## 5 Conditioning-Recogni…  -0.161   -0.296   -0.0265     1.96e-2   -0.590    0.268
## 6 Habituation-Recognit…  -0.127   -0.381    0.127      3.24e-1   -0.607    0.353

Type of reinforcement

If conditioning was used, was aversive or appetitive reinforcement used?

dat2 <- filter(dat, Type_reinforcement %in% c("Appetitive", "Aversive", "Not applicable"))
VCV_E2 <- impute_covariance_matrix(vi = dat2$lnRRV_E, cluster = dat2$Study_ID, r = 0.5)
mod_E3 <- rma.mv(yi = lnRR_Ea, V = VCV_E2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
                                                                                            ~1|ES_ID,
                                                                                            ~1|Strain),
                 test = "t",
                 data = dat2)

summary(mod_E3)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -7.8702   15.7405   27.7405   42.6723   28.7649   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0109  0.1046     30     no  Study_ID 
## sigma^2.2  0.0319  0.1787     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 764.5984, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.4138, p-val < .0001
## 
## Model Results:
## 
##                                   estimate      se    tval  df    pval    ci.lb 
## Type_reinforcementAppetitive        0.2175  0.0726  2.9942  89  0.0036   0.0732 
## Type_reinforcementAversive          0.2191  0.0410  5.3456  89  <.0001   0.1377 
## Type_reinforcementNot applicable    0.0812  0.0560  1.4504  89  0.1505  -0.0300 
##                                    ci.ub 
## Type_reinforcementAppetitive      0.3618   ** 
## Type_reinforcementAversive        0.3006  *** 
## Type_reinforcementNot applicable  0.1924      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_E3) 

##   R2_marginal R2_coditional 
##    0.07720103    1.00000000

Reinforcement_E <-orchard_plot(mod_E3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  xlim(-0.5, 2) + 
  scale_colour_manual(values = c("grey34","grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34", "grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
    geom_signif(comparisons = list(c("Aversive", "Not applicable")),
              annotations = "*") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Reinforcement_E

Fig. 4c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_E3 <- contrast_fun(data = dat2, response = lnRR_Ea, moderator = Type_reinforcement,
    VCV = VCV_E2)
res_table_mod_E3 <- get_estimate(model = mod_E3, contra = contra_mod_E3, moderator = Type_reinforcement)

res_table_mod_E3

## # A tibble: 6 x 7
##   Levels                 Estimate Lower_CI Upper_CI    P_value Lower_PI Upper_PI
##   <chr>                     <dbl>    <dbl>    <dbl>      <dbl>    <dbl>    <dbl>
## 1 Appetitive              0.217     0.0732   0.362     3.56e-3   -0.218    0.653
## 2 Aversive                0.219     0.138    0.301     6.89e-7   -0.200    0.638
## 3 Not applicable          0.0812   -0.0300   0.192     1.50e-1   -0.345    0.507
## 4 Appetitive-Aversive     0.00164  -0.163    0.167     9.84e-1   -0.442    0.445
## 5 Appetitive-Not applic… -0.136    -0.315    0.0425    1.33e-1   -0.585    0.312
## 6 Aversive-Not applicab… -0.138    -0.259   -0.0172    2.56e-2   -0.567    0.291

Age of enrichment

The age when individuals were exposed to environmental enrichment

dat3 <- filter(dat, Age_EE_exposure %in% c("Adult", "Adolescent"))
VCV_E3 <- impute_covariance_matrix(vi = dat3$lnRRV_E, cluster = dat3$Study_ID, r = 0.5)


mod_E4 <- rma.mv(yi = lnRR_Ea, V = VCV_E3, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
                                                                                    ~1|ES_ID,
                                                                                    ~1|Strain),
                 test = "t",
                 data = dat3)

summary(mod_E4) 

## 
## Multivariate Meta-Analysis Model (k = 88; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -6.9490   13.8980   23.8980   36.1698   24.6480   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0078  0.0885     29     no  Study_ID 
## sigma^2.2  0.0327  0.1808     88     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 86) = 782.6092, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 18.9060, p-val < .0001
## 
## Model Results:
## 
##                            estimate      se    tval  df    pval   ci.lb   ci.ub 
## Age_EE_exposureAdolescent    0.1799  0.0370  4.8553  86  <.0001  0.1062  0.2535 
## Age_EE_exposureAdult         0.2340  0.0620  3.7734  86  0.0003  0.1107  0.3573 
##  
## Age_EE_exposureAdolescent  *** 
## Age_EE_exposureAdult       *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_E4) 

##   R2_marginal R2_coditional 
##    0.01127347    1.00000000

Age_E<- orchard_plot(mod_E4, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  xlim(-0.5, 2) + 
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10))  

Age_E

Fig. 4d Orchard plot showing the group-wise means of the categorical variable ‘Age_EE_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_E4 <- contrast_fun(data = dat3, response = lnRR_Ea, moderator = Age_EE_exposure,
    VCV = VCV_E3)
res_table_mod_E4 <- get_estimate(model = mod_E4, contra = contra_mod_E4, moderator = Age_EE_exposure)

res_table_mod_E4

## # A tibble: 3 x 7
##   Levels           Estimate Lower_CI Upper_CI    P_value Lower_PI Upper_PI
##   <chr>               <dbl>    <dbl>    <dbl>      <dbl>    <dbl>    <dbl>
## 1 Adolescent         0.180    0.106     0.254 0.00000533   -0.227    0.587
## 2 Adult              0.234    0.111     0.357 0.000295     -0.185    0.653
## 3 Adolescent-Adult   0.0541  -0.0895    0.198 0.456        -0.371    0.479

Exercise enrichment

Did enrichment involve the addition of apparatus for voluntary exercise (e.g., a wheel or treadmill)?

mod_E5<- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~EE_exercise-1, random = list(~1|Study_ID,
                                                                               ~1|ES_ID,
                                                                               ~1|Strain),
                test = "t",
                data = dat)

summary(mod_E5)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -10.0303   20.0606   30.0606   42.5597   30.7749   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0096  0.0980     30     no  Study_ID 
## sigma^2.2  0.0345  0.1857     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 90) = 807.7169, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 16.1666, p-val < .0001
## 
## Model Results:
## 
##                         estimate      se    tval  df    pval   ci.lb   ci.ub 
## EE_exerciseExercise       0.1849  0.0407  4.5464  90  <.0001  0.1041  0.2657 
## EE_exerciseNo exercise    0.1900  0.0556  3.4151  90  0.0010  0.0795  0.3005 
##  
## EE_exerciseExercise     *** 
## EE_exerciseNo exercise  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_E5) 

##   R2_marginal R2_coditional 
##  0.0001360993  0.9999999987

Exercise_E <-orchard_plot(mod_E5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  xlim(-0.5, 2) + 
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Exercise_E

Fig. 4d Orchard plot showing the group-wise means of the categorical variable ‘EE_exercise’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_E5 <- contrast_fun(data = dat, response = lnRR_Ea, moderator = EE_exercise,
    VCV = VCV_E)
res_table_mod_E5 <- get_estimate(model = mod_E5, contra = contra_mod_E5, moderator = EE_exercise)

res_table_mod_E5

## # A tibble: 3 x 7
##   Levels               Estimate Lower_CI Upper_CI   P_value Lower_PI Upper_PI
##   <chr>                   <dbl>    <dbl>    <dbl>     <dbl>    <dbl>    <dbl>
## 1 Exercise              0.185     0.104     0.266 0.0000169   -0.240    0.610
## 2 No exercise           0.190     0.0795    0.301 0.000958    -0.241    0.621
## 3 Exercise-No exercise  0.00511  -0.132     0.142 0.941       -0.434    0.444

Social enrichment

Did enrichment involve more conspecifics (note that we did not include studies that only provided social enrichment but no other form of abiotic enrichment)?

dat4 <- filter(dat, EE_social %in% c("Social", "Non-social"))
VCV_E4 <- impute_covariance_matrix(vi = dat4$lnRRV_E, cluster = dat4$Study_ID, r = 0.5)
  
mod_E6<- rma.mv(yi = lnRR_Ea, V = VCV_E4, mod = ~EE_social-1, random = list(~1|Study_ID, 
                                                                             ~1|ES_ID,
                                                                             ~1|Strain),
                test = "t",data = dat4)

summary(mod_E6)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -9.3254   18.6507   28.6507   41.1498   29.3650   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0095  0.0974     30     no  Study_ID 
## sigma^2.2  0.0340  0.1844     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 90) = 801.5313, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 17.0110, p-val < .0001
## 
## Model Results:
## 
##                      estimate      se    tval  df    pval   ci.lb   ci.ub 
## EE_socialNon-social    0.1413  0.0498  2.8349  90  0.0057  0.0423  0.2403   ** 
## EE_socialSocial        0.2204  0.0432  5.0976  90  <.0001  0.1345  0.3063  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_E6) 

##   R2_marginal R2_coditional 
##    0.03377455    1.00000000

Social_E <-orchard_plot(mod_E6, mod = "EE_social", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  xlim(-0.5, 2) + 
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
       text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Social_E 

Fig. 4e Orchard plot showing the group-wise means of the categorical variable ‘EE_social’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_E6 <- contrast_fun(data = dat4, response = lnRR_Ea, moderator = EE_social,
    VCV = VCV_E4)
res_table_mod_E6 <- get_estimate(model = mod_E5, contra = contra_mod_E6, moderator = EE_social)

res_table_mod_E6

## # A tibble: 3 x 7
##   Levels                   Estimate Lower_CI Upper_CI  P_value Lower_PI Upper_PI
##   <chr>                       <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 EE_exerciseExercise        0.185    0.104     0.266  1.69e-5   -0.240    0.610
## 2 EE_exerciseNo exercise     0.190    0.0795    0.301  9.58e-4   -0.241    0.621
## 3 EE_exerciseExercise-EE_…   0.0791  -0.0520    0.210  2.34e-1   -0.355    0.514

Multi-moderator model

We ran a multi-moderator ‘full’ model of all moderators. We ran model selection based on AIC and selected the top set of models within delta AIC of 6 and calculated the importance of each moderator (i.e., Akaike weights) within this top model set.

# filter data so that all K < 5 are removed
dat_Efm <- dat %>%
    filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"), Type_reinforcement %in%
        c("Appetitive", "Aversive", "Not applicable"), EE_social %in% c("Social",
        "Non-social"), Age_EE_exposure %in% c("Adult", "Adolescent"))

# VCV_Efm <- impute_covariance_matrix(vi = dat_Efm$lnRRV_E, cluster =
# dat_Efm$Study_ID, r = 0.5)

mod_Efm <- rma.mv(yi = lnRR_Sa, V = VCV_Efm, mod = ~Type_assay - 1 + Learning_vs_memory +
    Type_reinforcement + EE_social + EE_exercise + Age_EE_exposure, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat_Efm)
# summary(mod_Efm) r2_ml(mod_Efm)

res_Efm <- dredge(mod_Efm, trace = 2)
saveRDS(res_Efm, file = here("Rdata", "res_Efm.rds"))
# also saving the full model and data
saveRDS(mod_Efm, file = here("Rdata", "mod_Efm.rds"))
saveRDS(dat_Efm, file = here("Rdata", "dat_Efm.rds"))

The akaike weights for the top set of models with AIC < 6

dat_Efm <- readRDS(file = here("Rdata", "dat_Efm.rds"))
mod_Efm <- readRDS(file = here("Rdata", "mod_Efm.rds"))
res_Efm <- readRDS(file = here("Rdata", "res_Efm.rds"))
res_Efm2 <- subset(res_Efm, delta <= 6, recalc.weights = FALSE)
importance(res_Efm2)

##                      Age_EE_exposure Type_assay EE_exercise EE_social
## Sum of weights:      0.63            0.36       0.16        0.14     
## N containing models:   11              10          5           5     
##                      Learning_vs_memory Type_reinforcement
## Sum of weights:      0.10               0.07              
## N containing models:    4                  4

Funnel plot

We produced funnel plots of effect sizes against the residuals of the full model to visually inspect for funnel asymmetry

Used to produce Fig. 7a

# funnel plot
Funnel_E <- funnel(mod_Efm, xlab = "lnRR", ylab = "Standard Error")

# Funnel_E

Fig.7a Funnel plot showing the standard error and residuals (lnRR) from the full model.

Egger’s regression and time-lag bias

We fitted effective sample size calculated as:

\[ \sqrt{\frac {1} { \tilde{N} }} = \sqrt{\frac {1} { N_\text{ES}} + \frac{1}{N_\text{EC}} + \frac {1}{N_\text{CS}} + \frac{1}{N_\text{CC}}}, \]

and year of publication as moderators in the full model to inferentially test for funnel asymmetry and a decrease in effect sizes with publication year

VCV_Efm <- impute_covariance_matrix(vi = dat_Efm$lnRRV_E, cluster = dat_Efm$Study_ID,
    r = 0.5)
dat_Efm$c_Year_published <- as.vector(scale(dat_Efm$Year_published, scale = F))

dat_Efm$sqrt_inv_e_n <- with(dat_Efm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))

PB_MR_E <- rma.mv(lnRR_Sa, VCV_Efm, mods = ~1 + sqrt_inv_e_n + Learning_vs_memory +
    Year_published + Type_assay + Type_reinforcement + EE_social + EE_exercise +
    Age_stress_exposure, random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML",
    test = "t", data = dat_Efm)

# PB_MR_E

estimates_PB_MR_E <- estimates.CI(PB_MR_E)
estimates_PB_MR_E

##                              estimate        mean        lower       upper
## 1                             intrcpt -9.48980951 -50.89864510 31.91902609
## 2                        sqrt_inv_e_n -0.06281057  -0.86857758  0.74295645
## 3            Learning_vs_memoryMemory -0.06858190  -0.17753633  0.04037252
## 4                      Year_published  0.00467568  -0.01591232  0.02526368
## 5               Type_assayHabituation -0.55467648  -0.99644702 -0.11290594
## 6               Type_assayRecognition -0.11960487  -0.45315282  0.21394309
## 7          Type_reinforcementAversive  0.15536918  -0.15356143  0.46429979
## 8    Type_reinforcementNot applicable  0.32741017  -0.11084706  0.76566739
## 9                     EE_socialSocial  0.04883161  -0.16642224  0.26408546
## 10             EE_exerciseNo exercise -0.03983196  -0.28097491  0.20131100
## 11           Age_stress_exposureAdult -0.16033122  -0.54043115  0.21976871
## 12 Age_stress_exposureEarly postnatal -0.05978986  -0.45195633  0.33237661
## 13        Age_stress_exposurePrenatal -0.13988902  -0.53943088  0.25965285

Leave-one-out analysis

We individually removed studies to determine if any singular studies hve a disproportionate effect on the meta-analytic mean and CI.

dat$Study_ID <- as.factor(dat$Study_ID)

LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$Study_ID))) {

    d <- dat %>%
        filter(Study_ID != levels(dat$Study_ID)[i])

    VCV_Eb <- impute_covariance_matrix(vi = d$lnRRV_E, cluster = d$Study_ID, r = 0.5)

    LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Ea, V = VCV_Eb, random = list(~1 |
        Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", data = dat[dat$Study_ID !=
        levels(dat$Study_ID)[i], ])
}

# TODO Shinichi to check

# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_E0) {
    df <- data.frame(est = mod_E0$b, lower = mod_E0$ci.lb, upper = mod_E0$ci.ub)
    return(df)
}



# using dplyr to form data frame
MA_CVR_E <- lapply(LeaveOneOut_effectsize, function(x) est.func(x)) %>%
    bind_rows %>%
    mutate(left_out = levels(dat$Study_ID))


saveRDS(MA_CVR_E, file = here("Rdata", "MA_CVR_E.rds"))

Used to produce Fig. S3

# telling ggplot to stop reordering factors
MA_CVR_E <- readRDS(file = here("Rdata", "MA_CVR_E.rds"))

MA_CVR_E$left_out <- as.factor(MA_CVR_E$left_out)
MA_CVR_E$left_out <- factor(MA_CVR_E$left_out, levels = MA_CVR_E$left_out)


# plotting
leaveoneout_E <- ggplot(MA_CVR_E) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +
    geom_hline(yintercept = mod_E0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
    geom_hline(yintercept = mod_E0$b, lty = 1, lwd = 0.75, colour = "black") + geom_hline(yintercept = mod_E0$ci.ub,
    lty = 3, lwd = 0.75, colour = "black") + geom_pointrange(aes(x = left_out, y = est,
    ymin = lower, ymax = upper)) + xlab("Study left out") + ylab("lnRR, 95% CI") +
    coord_flip() + theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
    theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))

leaveoneout_E

dat$Study_ID <- as.integer(dat$Study_ID)

Fig. S3 Leave-one-group-out analysis showing meta-analytic mean and 95% CI when each individual study is removed from the data set.

Using SMD

We re-ran the main effect of environmental enrichment meta-analytic model using standardised mean difference (SMD) instead of lnRR

mod_E0a <- rma.mv(yi = SMD_Ea, V = VCV_Ea, random = list(~1 | Study_ID, ~1 | ES_ID,
    ~1 | Strain), test = "t", data = dat)
summary(mod_E0a)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc 
## -111.6284   223.2568   231.2568   241.3003   231.7220   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.2512  0.5012     30     no  Study_ID 
## sigma^2.2  0.4247  0.6517     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Heterogeneity:
## Q(df = 91) = 673.4722, p-val < .0001
## 
## Model Results:
## 
## estimate      se    tval  df    pval   ci.lb   ci.ub 
##   0.7241  0.1295  5.5895  91  <.0001  0.4668  0.9814  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

i2_ml(mod_E0a)

##     I2_total  I2_Study_ID     I2_ES_ID    I2_Strain 
## 8.673320e-01 3.223202e-01 5.450118e-01 1.749563e-09

Risk of Bias

Including if the study was randomised as a moderator and conducting a pair-wise contrast

mod_randomised_E <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~ROB_randomisation - 1,
    random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)

summary(mod_randomised_E)
r2_ml(mod_randomised_E)

contra_mod_randomised_E <- contrast_fun(data = dat, response = lnRR_Ea, moderator = ROB_randomisation,
    VCV = VCV_E)
res_table_mod_randomised_E <- get_estimate(model = mod_randomised_E, contra = contra_mod_randomised_E,
    moderator = ROB_randomisation)

res_table_mod_randomised_E

## # A tibble: 3 x 7
##   Levels      Estimate Lower_CI Upper_CI    P_value Lower_PI Upper_PI
##   <chr>          <dbl>    <dbl>    <dbl>      <dbl>    <dbl>    <dbl>
## 1 Unclear       0.140    0.0488    0.232 0.00304      -0.276    0.556
## 2 Yes           0.227    0.139     0.315 0.00000158   -0.188    0.642
## 3 Unclear-Yes   0.0868  -0.0352    0.209 0.161        -0.337    0.511

Including if the study used blinding and conducting a pair-wise contrast

mod_blinding_E <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~ROB_blinding - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)

summary(mod_blinding_E)
r2_ml(mod_blinding_E)

contra_mod_blinding_E <- contrast_fun(data = dat, response = lnRR_Ea, moderator = ROB_blinding,
    VCV = VCV_E)
res_table_mod_blinding_E <- get_estimate(model = mod_blinding_E, contra = contra_mod_blinding_E,
    moderator = ROB_blinding)

res_table_mod_blinding_E

## # A tibble: 3 x 7
##   Levels      Estimate Lower_CI Upper_CI    P_value Lower_PI Upper_PI
##   <chr>          <dbl>    <dbl>    <dbl>      <dbl>    <dbl>    <dbl>
## 1 Unclear       0.198   0.122      0.274 0.00000131   -0.226    0.622
## 2 Yes           0.148   0.00811    0.288 0.0384       -0.292    0.588
## 3 Unclear-Yes  -0.0501 -0.208      0.108 0.530        -0.496    0.396

Learning vs Memory

Was the response learning or memory?

mod_S1 <-  rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
                                                                                        ~1|ES_ID,
                                                                                        ~1|Strain),
                  test = "t",
                  data = dat)

summary(mod_S1) 

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -14.5281   29.0562   39.0562   51.5553   39.7705   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0094  0.0970     30     no  Study_ID 
## sigma^2.2  0.0388  0.1969     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 90) = 946.3930, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 5.0145, p-val = 0.0086
## 
## Model Results:
## 
##                             estimate      se     tval  df    pval    ci.lb 
## Learning_vs_memoryLearning   -0.1211  0.0476  -2.5423  90  0.0127  -0.2157 
## Learning_vs_memoryMemory     -0.0974  0.0380  -2.5648  90  0.0120  -0.1728 
##                               ci.ub 
## Learning_vs_memoryLearning  -0.0265  * 
## Learning_vs_memoryMemory    -0.0219  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_S1) 

##   R2_marginal R2_coditional 
##   0.002776034   1.000000000

LvsM_S <- orchard_plot(mod_S1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
       text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

LvsM_S 

Fig. 5a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_S1 <- contrast_fun(data = dat, response = lnRR_Sa, moderator = Learning_vs_memory,
    VCV = VCV_S)
res_table_mod_S1 <- get_estimate(model = mod_S1, contra = contra_mod_S1, moderator = Learning_vs_memory)

res_table_mod_S1

## # A tibble: 3 x 7
##   Levels          Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>              <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Learning         -0.121   -0.216   -0.0265  0.0127   -0.567    0.325
## 2 Memory           -0.0974  -0.173   -0.0219  0.0120   -0.540    0.345
## 3 Learning-Memory   0.0237  -0.0780   0.125   0.644    -0.424    0.472

Type of assay

The broad category of the type of assay used to measure learning or memory

dat$Type_assay<-as.factor(dat$Type_assay)

VCV_S1 <- impute_covariance_matrix(vi = dat1$lnRRV_S, cluster = dat$Study_ID, r = 0.5)


mod_S2 <- rma.mv(yi = lnRR_Sa, V = VCV_S1, mod = ~Type_assay-1, random =   list(~1|Study_ID,
                                                                                    ~1|ES_ID,
                                                                                    ~1|Strain),
                 test = "t",
                 data = dat1)

summary(mod_S2)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -9.8028   19.6055   31.6055   46.5374   32.6299   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0161  0.1271     30     no  Study_ID 
## sigma^2.2  0.0279  0.1671     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 723.4973, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 6.7053, p-val = 0.0004
## 
## Model Results:
## 
##                         estimate      se     tval  df    pval    ci.lb    ci.ub 
## Type_assayConditioning   -0.0981  0.0375  -2.6192  89  0.0104  -0.1725  -0.0237 
## Type_assayHabituation    -0.4615  0.1126  -4.0969  89  <.0001  -0.6853  -0.2377 
## Type_assayRecognition    -0.0534  0.0645  -0.8287  89  0.4095  -0.1816   0.0747 
##  
## Type_assayConditioning    * 
## Type_assayHabituation   *** 
## Type_assayRecognition 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_S2) 

##   R2_marginal R2_coditional 
##     0.1853359     1.0000000

Learning_S <-orchard_plot(mod_S2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
    geom_signif(comparisons = list(c("Habituation", "Recognition")),
              annotations = "***") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Learning_S

Fig. 5b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_S2 <- contrast_fun(data = dat1, response = lnRR_Sa, moderator = Type_assay,
    VCV = VCV_S1)
res_table_mod_S2 <- get_estimate(model = mod_S2, contra = contra_mod_S2, moderator = Type_assay)

res_table_mod_S2

## # A tibble: 6 x 7
##   Levels                   Estimate Lower_CI Upper_CI  P_value Lower_PI Upper_PI
##   <chr>                       <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 Conditioning              -0.0981  -0.173   -0.0237  1.04e-2  -0.522    0.326 
## 2 Habituation               -0.461   -0.685   -0.238   9.21e-5  -0.935    0.0119
## 3 Recognition               -0.0534  -0.182    0.0747  4.10e-1  -0.490    0.383 
## 4 Conditioning-Habituation  -0.363   -0.584   -0.142   1.55e-3  -0.835    0.109 
## 5 Conditioning-Recognition   0.0447  -0.0882   0.177   5.06e-1  -0.393    0.482 
## 6 Habituation-Recognition    0.408    0.161    0.655   1.48e-3  -0.0768   0.893

Type of reinforcement

If conditioning was used, was aversive or appetitive reinforcement used?

VCV_S2 <- impute_covariance_matrix(vi = dat2$lnRRV_S, cluster = dat$Study_ID, r = 0.5)

mod_S3 <- rma.mv(yi = lnRR_Sa, V = VCV_S2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
                                                                                            ~1|ES_ID,
                                                                                            ~1|Strain),
                 test = "t",
                 data = dat2)

summary(mod_S3)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -13.8810   27.7621   39.7621   54.6939   40.7864   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0103  0.1016     30     no  Study_ID 
## sigma^2.2  0.0387  0.1966     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 920.8439, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.7942, p-val = 0.0130
## 
## Model Results:
## 
##                                   estimate      se     tval  df    pval 
## Type_reinforcementAppetitive       -0.1846  0.0749  -2.4649  89  0.0156 
## Type_reinforcementAversive         -0.0730  0.0427  -1.7081  89  0.0911 
## Type_reinforcementNot applicable   -0.1172  0.0590  -1.9851  89  0.0502 
##                                     ci.lb    ci.ub 
## Type_reinforcementAppetitive      -0.3334  -0.0358  * 
## Type_reinforcementAversive        -0.1579   0.0119  . 
## Type_reinforcementNot applicable  -0.2345   0.0001  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_S3) 

##   R2_marginal R2_coditional 
##     0.0366428     1.0000000

Reinforcement_S <-orchard_plot(mod_S3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
       text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Reinforcement_S

Fig. 5c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_S3 <- contrast_fun(data = dat2, response = lnRR_Sa, moderator = Type_reinforcement,
    VCV = VCV_S2)
res_table_mod_S3 <- get_estimate(model = mod_S3, contra = contra_mod_S3, moderator = Type_reinforcement)

res_table_mod_S3

## # A tibble: 6 x 7
##   Levels                    Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>                        <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Appetitive                 -0.185   -0.333  -3.58e-2  0.0156   -0.649    0.280
## 2 Aversive                   -0.0730  -0.158   1.19e-2  0.0911   -0.521    0.375
## 3 Not applicable             -0.117   -0.235   1.12e-4  0.0502   -0.572    0.338
## 4 Appetitive-Aversive         0.112   -0.0592  2.82e-1  0.198    -0.360    0.583
## 5 Appetitive-Not applicable   0.0674  -0.119   2.54e-1  0.474    -0.410    0.545
## 6 Aversive-Not applicable    -0.0442  -0.173   8.46e-2  0.497    -0.502    0.414

Age of stress

The age when individuals were exposed to stress

mod_S4 <-rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Age_stress_exposure-1, random = list(~1|Study_ID,
                                                                                       ~1|ES_ID,
                                                                                       ~1|Strain),
                test = "t",
                data = dat)
summary(mod_S4) 

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -12.4083   24.8166   38.8166   56.1579   40.2166   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0048  0.0694     30     no  Study_ID 
## sigma^2.2  0.0392  0.1979     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 88) = 881.1229, p-val < .0001
## 
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 4.3703, p-val = 0.0029
## 
## Model Results:
## 
##                                     estimate      se     tval  df    pval 
## Age_stress_exposureAdolescent         0.0074  0.1159   0.0641  88  0.9490 
## Age_stress_exposureAdult             -0.2279  0.0622  -3.6664  88  0.0004 
## Age_stress_exposureEarly postnatal   -0.0561  0.0435  -1.2891  88  0.2007 
## Age_stress_exposurePrenatal          -0.1145  0.0743  -1.5404  88  0.1271 
##                                       ci.lb    ci.ub 
## Age_stress_exposureAdolescent       -0.2228   0.2377      
## Age_stress_exposureAdult            -0.3514  -0.1044  *** 
## Age_stress_exposureEarly postnatal  -0.1425   0.0304      
## Age_stress_exposurePrenatal         -0.2621   0.0332      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_S4) 

##   R2_marginal R2_coditional 
##    0.09987307    1.00000000

Age_S <- orchard_plot(mod_S4, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
    geom_signif(comparisons = list(c("Adult", "Early postnatal")),
              annotations = "*") +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Age_S 

Fig. 5d Orchard plot showing the group-wise means of the categorical variable ‘Age_stress_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_S4 <- contrast_fun(data = dat, response = lnRR_Sa, moderator = Age_stress_exposure,
    VCV = VCV_S)
res_table_mod_S4 <- get_estimate(model = mod_S4, contra = contra_mod_S4, moderator = Age_stress_exposure)

res_table_mod_S4

## # A tibble: 10 x 7
##    Levels                   Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##    <chr>                       <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
##  1 Adolescent                0.00743  -0.223    0.238  9.49e-1   -0.469    0.484
##  2 Adult                    -0.228    -0.351   -0.104  4.20e-4   -0.663    0.207
##  3 Early postnatal          -0.0561   -0.142    0.0304 2.01e-1   -0.482    0.370
##  4 Prenatal                 -0.114    -0.262    0.0332 1.27e-1   -0.557    0.328
##  5 Adolescent-Adult         -0.235    -0.497    0.0260 7.70e-2   -0.727    0.257
##  6 Adolescent-Early postna… -0.0635   -0.309    0.182  6.09e-1   -0.547    0.421
##  7 Adolescent-Prenatal      -0.122    -0.395    0.152  3.78e-1   -0.620    0.377
##  8 Adult-Early postnatal     0.172     0.0211   0.323  2.60e-2   -0.271    0.615
##  9 Adult-Prenatal            0.113    -0.0791   0.306  2.45e-1   -0.346    0.573
## 10 Early postnatal-Prenatal -0.0584   -0.229    0.113  4.99e-1   -0.509    0.392

Type of stress

The type of stressor used

dat5 <- filter(dat, Type_stress_exposure %in% c("Restraint", "Noise", "MS", "Combination"))
VCV_S3 <- impute_covariance_matrix(vi = dat5$lnRRV_S, cluster = dat5$Study_ID, r = 0.5)

mod_S5 <- rma.mv(yi = lnRR_Sa, V = VCV_S3, mod = ~Type_stress_exposure-1, random = list(~1|Study_ID,
                                                                                         ~1|ES_ID,
                                                                                         ~1|Strain),
                 test = "t",
                 data = dat5)
summary(mod_S5) 

## 
## Multivariate Meta-Analysis Model (k = 85; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -11.4138   22.8276   36.8276   53.5887   38.3618   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0115  0.1071     25     no  Study_ID 
## sigma^2.2  0.0401  0.2002     85     no     ES_ID 
## sigma^2.3  0.0000  0.0000      4     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 81) = 897.4553, p-val < .0001
## 
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 81) = 2.8767, p-val = 0.0278
## 
## Model Results:
## 
##                                  estimate      se     tval  df    pval    ci.lb 
## Type_stress_exposureCombination   -0.0500  0.0892  -0.5605  81  0.5767  -0.2274 
## Type_stress_exposureMS            -0.0539  0.0560  -0.9630  81  0.3384  -0.1653 
## Type_stress_exposureNoise         -0.1203  0.1036  -1.1608  81  0.2491  -0.3265 
## Type_stress_exposureRestraint     -0.2101  0.0704  -2.9863  81  0.0037  -0.3501 
##                                    ci.ub 
## Type_stress_exposureCombination   0.1274     
## Type_stress_exposureMS            0.0575     
## Type_stress_exposureNoise         0.0859     
## Type_stress_exposureRestraint    -0.0701  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_S5)

##   R2_marginal R2_coditional 
##    0.07029554    1.00000000

Stressor<- orchard_plot(mod_S5, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Stressor

Fig. 5e Orchard plot showing the group-wise means of the categorical variable ‘Type_stress_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_S5 <- contrast_fun(data = dat5, response = lnRR_Sa, moderator = Type_stress_exposure,
    VCV = VCV_S3)
res_table_mod_S5 <- get_estimate(model = mod_S5, contra = contra_mod_S5, moderator = Type_stress_exposure)

res_table_mod_S5

## # A tibble: 10 x 7
##    Levels                Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##    <chr>                    <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
##  1 Combination           -0.0500    -0.227   0.127  0.577     -0.535    0.435
##  2 MS                    -0.0539    -0.165   0.0575 0.338     -0.519    0.411
##  3 Noise                 -0.120     -0.326   0.0859 0.249     -0.617    0.376
##  4 Restraint             -0.210     -0.350  -0.0701 0.00373   -0.683    0.263
##  5 Combination-MS        -0.00393   -0.213   0.206  0.970     -0.502    0.494
##  6 Combination-Noise     -0.0703    -0.342   0.202  0.608     -0.598    0.457
##  7 Combination-Restraint -0.160     -0.386   0.0658 0.162     -0.665    0.345
##  8 MS-Noise              -0.0664    -0.301   0.168  0.575     -0.575    0.443
##  9 MS-Restraint          -0.156     -0.335   0.0227 0.0861    -0.642    0.330
## 10 Noise-Restraint       -0.0898    -0.339   0.159  0.475     -0.606    0.426

Stess duration

Was the stress acute or chronic?

dat6 <- filter(dat, Stress_duration %in% c("Chronic", "Acute"))
VCV_S4 <- impute_covariance_matrix(vi = dat6$lnRRV_S, cluster = dat6$Study_ID, r = 0.5)

mod_S6 <-rma.mv(yi = lnRR_Sa, V = VCV_S4, mod = ~Stress_duration-1, random = list(~1|Study_ID,
                                                                                   ~1|ES_ID,
                                                                                   ~1|Strain),
                test = "t",
                data = dat6)
summary(mod_S6) 

## 
## Multivariate Meta-Analysis Model (k = 89; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -13.7898   27.5796   37.5796   49.9092   38.3204   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0062  0.0786     29     no  Study_ID 
## sigma^2.2  0.0391  0.1978     89     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 87) = 915.9393, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 87) = 6.6661, p-val = 0.0020
## 
## Model Results:
## 
##                         estimate      se     tval  df    pval    ci.lb    ci.ub 
## Stress_durationAcute      0.0135  0.0659   0.2045  87  0.8384  -0.1176   0.1445 
## Stress_durationChronic   -0.1360  0.0373  -3.6456  87  0.0005  -0.2101  -0.0618 
##  
## Stress_durationAcute 
## Stress_durationChronic  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_S6) 

##   R2_marginal R2_coditional 
##    0.08245896    1.00000000

Duration_S <- orchard_plot(mod_S6, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
    geom_signif(comparisons = list(c("Chronic", "Acute")),
              annotations = "*") +
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 15), # change font sizes
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10)) 

Duration_S 

Fig. 5f Orchard plot showing the group-wise means of the categorical variable ‘Stress_duration’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_S6 <- contrast_fun(data = dat6, response = lnRR_Sa, moderator = Stress_duration,
    VCV = VCV_S4)
res_table_mod_S6 <- get_estimate(model = mod_S6, contra = contra_mod_S6, moderator = Stress_duration)

res_table_mod_S6

## # A tibble: 3 x 7
##   Levels        Estimate Lower_CI Upper_CI  P_value Lower_PI Upper_PI
##   <chr>            <dbl>    <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 Acute           0.0135   -0.118  0.145   0.838      -0.429    0.456
## 2 Chronic        -0.136    -0.210 -0.0618  0.000453   -0.565    0.294
## 3 Acute-Chronic  -0.149    -0.300  0.00111 0.0517     -0.599    0.300

Multi-moderator model

We ran a multi-moderator ‘full’ model of all moderators. We ran model selection based on AIC and selected the top set of models within delta AIC of 6 and calculated the importance of each moderator (i.e., Akaike weights) within this top model set.

# selecting moderator levels with k >=5
dat_Sfm <- dat %>%
    filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"), Type_reinforcement %in%
        c("Appetitive", "Aversive", "Not applicable"), Type_stress_exposure %in%
        c("Restraint", "Noise", "MS", "Combination"), Stress_duration %in% c("Chronic",
        "Acute"))

# VCV_Sfm <- impute_covariance_matrix(vi = dat_Sfm$lnRRV_E, cluster =
# dat_Sfm$Study_ID, r = 0.5)

mod_Sfm <- rma.mv(yi = lnRR_Sa, V = VCV_Sfm, mod = ~Type_assay - 1 + Learning_vs_memory +
    Type_reinforcement + Type_stress_exposure + Age_stress_exposure + Stress_duration,
    random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat_Sfm)
# summary(mod_Sfm) r2_ml(mod_Sfm)

res_Sfm <- dredge(mod_Sfm, trace = 2)
saveRDS(res_Sfm, file = here("Rdata", "res_Sfm.rds"))
# also saving the full model and data
saveRDS(mod_Sfm, file = here("Rdata", "mod_Sfm.rds"))
saveRDS(dat_Sfm, file = here("Rdata", "dat_Sfm.rds"))

The akaike weights for the top set of models with AIC < 6

dat_Sfm <- readRDS(file = here("Rdata", "dat_Sfm.rds"))
mod_Sfm <- readRDS(file = here("Rdata", "mod_Sfm.rds"))
res_Sfm <- readRDS(file = here("Rdata", "res_Sfm.rds"))
res_Sfm2 <- subset(res_Sfm, delta <= 6, recalc.weights = FALSE)
importance(res_Sfm2)

##                      Type_assay Stress_duration Type_stress_exposure
## Sum of weights:      0.91       0.88            0.22                
## N containing models:    7          6               2                
##                      Learning_vs_memory Age_stress_exposure Type_reinforcement
## Sum of weights:      0.16               0.04                0.04              
## N containing models:    2                  1                   1

Funnel plot

We produced funnel plots of effect sizes against the residuals of the full model to visually inspect for funnel asymmetry.

Used to produce Fig. 7b

# funnel plot
Funnel_S <- funnel(mod_Sfm, xlab = "lnRR", ylab = "Standard Error")

# Funnel_S

Fig.7b Funnel plot showing the standard error and residuals (lnRR) from the full model.

Egger’s regression and time-lag bias

We fitted effective sample size calculated as:

\[ \sqrt{\frac {1} { \tilde{N} }} = \sqrt{\frac {1} { N_\text{ES}} + \frac{1}{N_\text{EC}} + \frac {1}{N_\text{CS}} + \frac{1}{N_\text{CC}}}, \]

and year of publication as moderators in the full model to inferentially test for funnel asymmetry and a decrease in effect sizes with publication year.

# calculating inv effective sample size for use in full meta-regression
VCV_Sfm <- impute_covariance_matrix(vi = dat_Sfm$lnRRV_E, cluster = dat_Sfm$Study_ID,
    r = 0.5)

dat_Sfm$sqrt_inv_e_n <- with(dat_Sfm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))

# time lag bias and eggers regression
dat_Sfm$c_Year_published <- as.vector(scale(dat_Sfm$Year_published, scale = F))

PB_MR_S <- rma.mv(lnRR_Sa, VCV_Sfm, mods = ~1 + sqrt_inv_e_n + c_Year_published +
    Type_assay + Learning_vs_memory + Type_reinforcement + Type_stress_exposure +
    Age_stress_exposure + Stress_duration, random = list(~1 | Study_ID, ~1 | ES_ID,
    ~1 | Strain), method = "REML", test = "t", data = dat_Sfm, control = list(optimizer = "optim",
    optmethod = "Nelder-Mead"))

# PB_MR_S

estimates_PB_MR_S <- estimates.CI(PB_MR_S)
estimates_PB_MR_S

##                              estimate         mean       lower       upper
## 1                             intrcpt  0.208757092 -0.47118467  0.88869885
## 2                        sqrt_inv_e_n  0.109877719 -0.90636737  1.12612281
## 3                    c_Year_published -0.006425821 -0.02804904  0.01519740
## 4               Type_assayHabituation -0.663428076 -1.09349872 -0.23335743
## 5               Type_assayRecognition -0.069789046 -0.41156145  0.27198335
## 6            Learning_vs_memoryMemory -0.076712961 -0.19562432  0.04219839
## 7          Type_reinforcementAversive  0.039563987 -0.13761334  0.21674132
## 8    Type_reinforcementNot applicable  0.189481532 -0.16955283  0.54851589
## 9              Type_stress_exposureMS  0.233030567 -0.40478886  0.87084999
## 10          Type_stress_exposureNoise -0.020700600 -0.76667312  0.72527192
## 11      Type_stress_exposureRestraint -0.067490308 -0.43366600  0.29868539
## 12           Age_stress_exposureAdult -0.129875123 -0.51040393  0.25065368
## 13 Age_stress_exposureEarly postnatal -0.193461086 -0.89857848  0.51165631
## 14        Age_stress_exposurePrenatal -0.146406300 -0.57377820  0.28096560
## 15             Stress_durationChronic -0.395711964 -0.65566557 -0.13575836

Leave-one-out sensitivity analysis

We individually removed studies to determine if any singular studies hve a disproportionate effect on the meta-analytic mean and CI.

dat$Study_ID <- as.factor(dat$Study_ID)

LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$Study_ID))) {
    d <- dat %>%
        filter(Study_ID != levels(dat$Study_ID)[i])

    VCV_Sb <- impute_covariance_matrix(vi = d$lnRRV_E, cluster = d$Study_ID, r = 0.5)

    LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Sa, V = VCV_Sb, random = list(~1 |
        Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", data = dat[dat$Study_ID !=
        levels(dat$Study_ID)[i], ])
}


# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_S0) {
    df <- data.frame(est = mod_S0$b, lower = mod_S0$ci.lb, upper = mod_S0$ci.ub)
    return(df)
}


# using dplyr to form data frame
MA_CVR_S <- lapply(LeaveOneOut_effectsize, function(x) est.func(x)) %>%
    bind_rows %>%
    mutate(left_out = levels(dat$Study_ID))

saveRDS(MA_CVR_S, file = here("Rdata", "MA_CVR_S.rds"))

Used to produce Fig. S4

MA_CVR_S <- readRDS(file = here("Rdata", "MA_CVR_S.rds"))

# telling ggplot to stop reordering factors
MA_CVR_S$left_out <- as.factor(MA_CVR_S$left_out)
MA_CVR_S$left_out <- factor(MA_CVR_S$left_out, levels = MA_CVR_S$left_out)

# plotting
leaveoneout_S <- ggplot(MA_CVR_S) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +
    geom_hline(yintercept = mod_S0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
    geom_hline(yintercept = mod_S0$b, lty = 1, lwd = 0.75, colour = "black") + geom_hline(yintercept = mod_S0$ci.ub,
    lty = 3, lwd = 0.75, colour = "black") + geom_pointrange(aes(x = left_out, y = est,
    ymin = lower, ymax = upper)) + xlab("Study left out") + ylab("lnRR, 95% CI") +
    coord_flip() + theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
    theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))

leaveoneout_S

dat$Study_ID <- as.integer(dat$Study_ID)

Fig. S4 Leave-one-group-out analysis showing meta-analytic mean and 95% CI when each individual study is removed from the data set.

Using SMD

We re-ran the main effect of stress meta-analytic model using standardised mean difference (SMD) instead of lnRR

mod_S0a <- rma.mv(yi = SMD_Sa, V = VCV_Sa, random = list(~1 | Study_ID, ~1 | ES_ID,
    ~1 | Strain), test = "t", data = dat)
summary(mod_S0a)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc 
## -120.2817   240.5634   248.5634   258.6068   249.0285   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.3874  0.6224     30     no  Study_ID 
## sigma^2.2  0.4969  0.7049     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Heterogeneity:
## Q(df = 91) = 818.5592, p-val < .0001
## 
## Model Results:
## 
## estimate      se     tval  df    pval    ci.lb    ci.ub 
##  -0.4256  0.1498  -2.8403  91  0.0056  -0.7232  -0.1279  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

i2_ml(mod_S0a)

##     I2_total  I2_Study_ID     I2_ES_ID    I2_Strain 
## 8.959540e-01 3.924790e-01 5.034750e-01 1.438454e-09

Risk of Bias

Including if the study was randomised as a moderator and conducting a pair-wise contrast

mod_randomised_S <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~ROB_randomisation - 1,
    random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)

summary(mod_randomised_S)
r2_ml(mod_randomised_S)

contra_mod_randomised_S <- contrast_fun(data = dat, response = lnRR_Sa, moderator = ROB_randomisation,
    VCV = VCV_S)
res_table_mod_randomised_S <- get_estimate(model = mod_randomised_S, contra = contra_mod_randomised_S,
    moderator = ROB_randomisation)

res_table_mod_randomised_S

## # A tibble: 3 x 7
##   Levels      Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>          <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Unclear      -0.0672   -0.164   0.0291 0.169     -0.509    0.374
## 2 Yes          -0.139    -0.230  -0.0480 0.00317   -0.580    0.301
## 3 Unclear-Yes  -0.0719   -0.204   0.0607 0.284     -0.523    0.379

Including if the study used blinding and conducting a pair-wise contrast

mod_blinding_S <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~ROB_blinding - 1, random = list(~1 |
    Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat)

summary(mod_blinding_S)
r2_ml(mod_blinding_S)

contra_mod_blinding_S <- contrast_fun(data = dat, response = lnRR_Sa, moderator = ROB_blinding,
    VCV = VCV_S)
res_table_mod_blinding_S <- get_estimate(model = mod_blinding_S, contra = contra_mod_blinding_S,
    moderator = ROB_blinding)

res_table_mod_blinding_S

## # A tibble: 3 x 7
##   Levels      Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>          <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Unclear      -0.119    -0.196  -0.0419 0.00288   -0.564    0.326
## 2 Yes          -0.0567   -0.202   0.0890 0.442     -0.518    0.405
## 3 Unclear-Yes   0.0625   -0.102   0.227  0.454     -0.406    0.531

Learning vs Memory

Was the response learning or memory?

mod_ES1 <-  rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID, 
                                                                                           ~1|ES_ID,
                                                                                           ~1|Strain),
                   test = "t",
                   data = dat)

summary(mod_ES1) 

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -40.4769   80.9539   90.9539  103.4529   91.6682   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0292  0.1708     30     no  Study_ID 
## sigma^2.2  0.0232  0.1524     92     no     ES_ID 
## sigma^2.3  0.0034  0.0582      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 90) = 299.1854, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 2.9219, p-val = 0.0590
## 
## Model Results:
## 
##                             estimate      se    tval  df    pval    ci.lb 
## Learning_vs_memoryLearning    0.1744  0.0722  2.4166  90  0.0177   0.0310 
## Learning_vs_memoryMemory      0.1057  0.0619  1.7065  90  0.0914  -0.0174 
##                              ci.ub 
## Learning_vs_memoryLearning  0.3178  * 
## Learning_vs_memoryMemory    0.2288  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_ES1) 

##   R2_marginal R2_coditional 
##     0.0197648     0.9405080

LvsM_ES <- orchard_plot(mod_ES1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 3, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        axis.text.x = element_text(size = 10), # change font sizes
        axis.text.y = element_text(size = 10),
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 7)) 

LvsM_ES 

Fig. 6a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_ES1 <- contrast_fun(data = dat, response = lnRR_ESa, moderator = Learning_vs_memory,
    VCV = VCV_ES)
res_table_mod_ES1 <- get_estimate(model = mod_ES1, contra = contra_mod_ES1, moderator = Learning_vs_memory)

res_table_mod_ES1

## # A tibble: 3 x 7
##   Levels          Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>              <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Learning          0.174    0.0310   0.318   0.0177   -0.316    0.665
## 2 Memory            0.106   -0.0174   0.229   0.0914   -0.379    0.591
## 3 Learning-Memory  -0.0687  -0.176    0.0382  0.205    -0.550    0.413

Type of assay

The broad category of the type of assay used to measure learning or memory

VCV_ES1 <- impute_covariance_matrix(vi = dat1$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)

mod_ES2 <- rma.mv(yi = lnRR_ESa, V = VCV_ES1, mod = ~Type_assay-1, random = list(~1|Study_ID,
                                                                                    ~1|ES_ID,
                                                                                    ~1|Strain),
                  test = "t",
                  data = dat1)

summary(mod_ES2)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -39.0874   78.1748   90.1748  105.1066   91.1992   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0370  0.1924     30     no  Study_ID 
## sigma^2.2  0.0192  0.1386     92     no     ES_ID 
## sigma^2.3  0.0018  0.0422      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 293.9385, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.1062, p-val = 0.0305
## 
## Model Results:
## 
##                         estimate      se     tval  df    pval    ci.lb   ci.ub 
## Type_assayConditioning    0.1525  0.0589   2.5877  89  0.0113   0.0354  0.2696 
## Type_assayHabituation     0.1990  0.1415   1.4070  89  0.1629  -0.0820  0.4801 
## Type_assayRecognition    -0.0048  0.0800  -0.0606  89  0.9518  -0.1637  0.1541 
##  
## Type_assayConditioning  * 
## Type_assayHabituation 
## Type_assayRecognition 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_ES2) 

##   R2_marginal R2_coditional 
##    0.05775809    0.97105550

Learning_ES <- orchard_plot(mod_ES2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34")) +
  geom_point(aes(fill = name),  size = 3, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
    geom_signif(comparisons = list(c("Conditioning", "Recognition")),
              annotations = "*") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        axis.text.x = element_text(size = 10), # change font sizes
        axis.text.y = element_text(size = 10),
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 7)) 

Learning_ES

Fig. 6b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_ES2 <- contrast_fun(data = dat1, response = lnRR_ESa, moderator = Type_assay,
    VCV = VCV_ES1)
res_table_mod_ES2 <- get_estimate(model = mod_ES2, contra = contra_mod_ES2, moderator = Type_assay)

res_table_mod_ES2

## # A tibble: 6 x 7
##   Levels                   Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>                       <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Conditioning              0.153     0.0354   0.270   0.0113   -0.340    0.645
## 2 Habituation               0.199    -0.0820   0.480   0.163    -0.356    0.754
## 3 Recognition              -0.00485  -0.164    0.154   0.952    -0.509    0.499
## 4 Conditioning-Habituation  0.0465   -0.219    0.312   0.728    -0.501    0.594
## 5 Conditioning-Recognition -0.157    -0.304   -0.0105  0.0360   -0.658    0.343
## 6 Habituation-Recognition  -0.204    -0.488    0.0801  0.157    -0.760    0.353

Type of reinforcement

If conditioning was used, was aversive or appetitive reinforcement used?

VCV_ES2 <- impute_covariance_matrix(vi = dat2$lnRRV_ES, cluster = dat2$Study_ID, r = 0.5)

mod_ES3 <- rma.mv(yi = lnRR_ESa, V = VCV_ES2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
                                                                                               ~1|ES_ID,
                                                                                               ~1|Strain),
                  test = "t",
                  data = dat2)

summary(mod_ES3)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -39.0604   78.1208   90.1208  105.0526   91.1452   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0382  0.1954     30     no  Study_ID 
## sigma^2.2  0.0189  0.1377     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 293.4724, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.2547, p-val = 0.0254
## 
## Model Results:
## 
##                                   estimate      se    tval  df    pval    ci.lb 
## Type_reinforcementAppetitive        0.1007  0.1075  0.9366  89  0.3515  -0.1129 
## Type_reinforcementAversive          0.1573  0.0569  2.7618  89  0.0070   0.0441 
## Type_reinforcementNot applicable    0.0147  0.0702  0.2101  89  0.8341  -0.1247 
##                                    ci.ub 
## Type_reinforcementAppetitive      0.3143     
## Type_reinforcementAversive        0.2704  ** 
## Type_reinforcementNot applicable  0.1541     
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_ES3) 

##   R2_marginal R2_coditional 
##     0.0586952     1.0000000

Reinforcement_ES <- orchard_plot(mod_ES3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34","grey34")) +
  geom_point(aes(fill = name),  size = 3, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
    geom_signif(comparisons = list(c("Aversive", "Not applicable")),
              annotations = "*") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        axis.text.x = element_text(size = 10), # change font sizes
        axis.text.y = element_text(size = 10),
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 7)) 

Reinforcement_ES 

Fig. 6c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_ES3 <- contrast_fun(data = dat2, response = lnRR_ESa, moderator = Type_reinforcement,
    VCV = VCV_ES2)
res_table_mod_ES3 <- get_estimate(model = mod_ES3, contra = contra_mod_ES3, moderator = Type_reinforcement)

res_table_mod_ES3

## # A tibble: 6 x 7
##   Levels                    Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>                        <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Appetitive                  0.101   -0.113   0.314   0.351     -0.420    0.621
## 2 Aversive                    0.157    0.0441  0.270   0.00698   -0.331    0.645
## 3 Not applicable              0.0147  -0.125   0.154   0.834     -0.480    0.510
## 4 Appetitive-Aversive         0.0566  -0.182   0.295   0.639     -0.475    0.588
## 5 Appetitive-Not applicable  -0.0859  -0.333   0.161   0.492     -0.621    0.450
## 6 Aversive-Not applicable    -0.143   -0.278  -0.00744 0.0389    -0.636    0.351

Age of enrichment

The age when individuals were exposed to environmental enrichment

VCV_ES3 <- impute_covariance_matrix(vi = dat3$lnRRV_ES, cluster = dat3$Study_ID, r = 0.5)

mod_ES4 <- rma.mv(yi = lnRR_ESa, V = VCV_ES3, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
                                                                                    ~1|ES_ID,
                                                                                    ~1|Strain),
                 test = "t",
                 data = dat3)

summary(mod_ES4) 

## 
## Multivariate Meta-Analysis Model (k = 88; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -36.6407   73.2813   83.2813   95.5531   84.0313   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0332  0.1822     29     no  Study_ID 
## sigma^2.2  0.0207  0.1439     88     no     ES_ID 
## sigma^2.3  0.0007  0.0265      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 86) = 288.6493, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 3.1308, p-val = 0.0487
## 
## Model Results:
## 
##                            estimate      se    tval  df    pval    ci.lb 
## Age_EE_exposureAdolescent    0.1284  0.0583  2.2006  86  0.0304   0.0124 
## Age_EE_exposureAdult         0.1247  0.0921  1.3538  86  0.1793  -0.0584 
##                             ci.ub 
## Age_EE_exposureAdolescent  0.2444  * 
## Age_EE_exposureAdult       0.3077    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_ES4) 

##   R2_marginal R2_coditional 
##  0.0000400928  0.9871095822

Age_enrichment_ES <- orchard_plot(mod_ES4, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 3, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
       axis.text.x = element_text(size = 10), # change font sizes
        axis.text.y = element_text(size = 10),
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 7)) 

Age_enrichment_ES

Fig. 6d Orchard plot showing the group-wise means of the categorical variable ‘Age_EE_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_ES4 <- contrast_fun(data = dat3, response = lnRR_ESa, moderator = Age_EE_exposure,
    VCV = VCV_ES3)
res_table_mod_ES4 <- get_estimate(model = mod_ES4, contra = contra_mod_ES4, moderator = Age_EE_exposure)

res_table_mod_ES4

## # A tibble: 3 x 7
##   Levels           Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>               <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Adolescent        0.128     0.0124    0.244  0.0304   -0.350    0.607
## 2 Adult             0.125    -0.0584    0.308  0.179    -0.375    0.624
## 3 Adolescent-Adult -0.00373  -0.213     0.205  0.972    -0.513    0.506

Exercise enrichment

Did enrichment involve the addition of apparatus for voluntary exercise (e.g., a wheel or treadmill)?

mod_ES5<- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~EE_exercise-1, random = list(~1|Study_ID, 
    ~1|ES_ID,
    ~1|Strain),
     test = "t",
     data = dat)

summary(mod_ES5)

## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -40.4968   80.9935   90.9935  103.4926   91.7078   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0325  0.1803     30     no  Study_ID 
## sigma^2.2  0.0230  0.1517     92     no     ES_ID 
## sigma^2.3  0.0065  0.0805      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 90) = 297.3270, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 1.8401, p-val = 0.1647
## 
## Model Results:
## 
##                         estimate      se    tval  df    pval    ci.lb   ci.ub 
## EE_exerciseExercise       0.1051  0.0780  1.3474  90  0.1812  -0.0499  0.2602 
## EE_exerciseNo exercise    0.1687  0.0967  1.7448  90  0.0844  -0.0234  0.3609 
##  
## EE_exerciseExercise 
## EE_exerciseNo exercise  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

r2_ml(mod_ES5) 

##   R2_marginal R2_coditional 
##    0.01474459    0.89712469

Exercise_ES <- orchard_plot(mod_ES5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  scale_colour_manual(values = c("grey34","grey34")) +
  scale_fill_manual(values=c("grey34","grey34")) +
  geom_point(aes(fill = name),  size = 3, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
       axis.text.x = element_text(size = 10), # change font sizes
        axis.text.y = element_text(size = 10),
        legend.title = element_text(size = 7),
        legend.text = element_text(size = 7)) 

Exercise_ES

Fig. 6d Orchard plot showing the group-wise means of the categorical variable ‘EE_exercise’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.

contra_mod_ES5 <- contrast_fun(data = dat, response = lnRR_ESa, moderator = EE_exercise,
    VCV = VCV_ES)
res_table_mod_ES5 <- get_estimate(model = mod_ES5, contra = contra_mod_ES5, moderator = EE_exercise)

res_table_mod_ES5

## # A tibble: 3 x 7
##   Levels               Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
##   <chr>                   <dbl>    <dbl>    <dbl>   <dbl>    <dbl>    <dbl>
## 1 Exercise               0.105   -0.0499    0.260  0.181    -0.413    0.624
## 2 No exercise            0.169   -0.0234    0.361  0.0844   -0.362    0.699
## 3 Exercise-No exercise   0.0636  -0.138     0.265  0.532    -0.470    0.598

Social enrichment

Did enrichment involve more conspecifics (note that we did not include studies that only provided social enrichment but no other form of abiotic enrichment)?

VCV_ES4 <- impute_covariance_matrix(vi = dat4$lnRRV_ES, cluster = dat4$Study_ID, r = 0.5)
mod_ES6<- rma.mv(yi = lnRR_ESa, V = VCV_ES4, mod = ~EE_social-1, random = list(~1|Study_ID,
                                                                                ~1|ES_ID,
                                                                                ~1|Strain),
                 test = "t",
                 data = dat4)

summary(mod_ES6)